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We analyze recent experiments on the dilute rare-earth compound LiHozYi-a^ in the context 
of an effective Ising dipolar model. Using a Monte Carlo method we calculate the low-temperature 
behavior of the specific heat and linear susceptibility, and compare our results to measurements. In 
our model the susceptibility follows a Curie- Weiss law at high temperature, \ ~ ^-/iT — Tew), with a 
Curie- Weiss temperature that scales with dilution, Tew ~ x, consistent with early experiments. We 
also find that the peak in the specific heat scales linearly with dilution, Cmax (T) ~ x, in disagreement 
with recent experiments. Experimental studies do not reach a consensus on the functional form 
of these quantities, and in particular we do not see reported scalings of the form \ ~ T~~ ' 75 and 
X ~ exp (— T/To). Furthermore we calculate the ground state magnetization as a function of dilution, 
and re-examine the phase diagram around the critical dilution x c = 0.24 ± 0.03. We find that the 
spin glass susceptibility for the Ising model does not diverge below x c , while recent experiments give 
strong evidence for a stable spin-glass phase in L1H00.167Y0.833F4. 

PACS numbers: 75.10.Hk,75.50.Lk,75.40.Mg 



I. INTRODUCTION 

The rare-earth compound LiHoF4 is used as a model 
magnet to investigate diverse magnetic phenomena such 
as quantum phase transitions, 1 spin-glass behavior— and 
quantum annealing 3 -. The magnetic behavior arises from 
the Ho 3+ ions which have tightly bound 4/ electrons. 
This causes the exchange interaction to be weak, and the 
inter- ion interactions are predominantly dipolar. The lo- 
cal crystal field causes a strong anisotropy, and the in- 
teraction is Ising-like. To a first approximation LiHoF4 
is therefore believed to be good realization of a dipolar 
Ising modelj^ 

i^j l 3 i,nn 

where we have used a dipolar coupling constant J = 
0.214 K and a nearest-neighbor (nn) exchange coupling 
Jcx = 0.12 K£ The interspin distance is ry, with a com- 
ponent Zij along the Ising axis. The magnetic Ho 3+ 
ions sit on a tetragonal lattice with four ions per unit 
cell. To study quantum criticality a transverse magnetic 
field can be applied, and in order to study the effects of 
disorder the magnetic Ho 3+ ions can be substituted by 
nonmagnetic Y 3+ ions, resulting in LiHo a; Yi_ a; F4. Dur- 
ing the last three decades LHI0F4 has been extensively 
studied and used as a textbook example of a quantum 
magnet£i However, in the case of substantial dilution 
measurements have reported a variety of functional forms 
for basic thermodynamic quantities, such as the static 
susceptibility and the specific heat. 

The earliest data we find for the static susceptibil- 
ity report a high-temperature Curie- Weiss scaling x ~ 
1/(T — T cw ) with Curie- Weiss temperatures Tew = 0.05 
and 0.16 for dilution x = 0.045 and 0.167 respectively.™ 
In a later work the susceptibility is found to diverge with 
a different power law, x ~ 7 1-0 - 75 (x = 0.045) and 



in a recent study the exponential low-temperature form 
X = cxp(— T/Tq) is reported^ The specific heat has also 
been measured by several different groups, and in an ear- 
lier study^ of the specific heat a peak was found at about 
T = 0.3 K for x = 0.045, while there was only a much 
broader maximum below T = 0.2 K for x=0.167. 2 A later 
study by the same group finds peaks at T = 0.1 K, as well 
as at T = 0.3 K, for x = 0.045. 8 A more recent studyiS 
displays a dilution independent maximum in the specific 
heat at about T = 0.1 K for x = 0.018, 0.045 and 0.08. 

Finally, the nature of the glassy phase at low tem- 
peratures has also been the topic of several experimen- 
tal studies. Earlier work found a spin-liquid (anti-glass) 
phase at extreme dilution (x=0.045), followed by a stable 
spin-glass phase at dilution x=0.167, and finally a mag- 
netic phase at x=0.3i£ More recent experiments did not 
detect a spin glass transition^ but this may have been 
due to the use of large magnetic fields p^i Recent numer- 
ical work on dilute dipoles on a small cubic lattice fails 
to find a spin glass transition, 1 - and so does a recent 
numerical study of the above model for LiHo :E Yi_ a ;F4i^ 

In this study we confine ourselves to the case of no 
external magnetic field, but it is interesting to note that 
quantum Monte Carlo studies of the above non-diluted 
model including an applied transverse field^ 3 - do not 
reach good agreement with the experimental phase di- 
agram, even for small transverse fields. 

Much of the recent theoretical work on LiHo x Yi_ x F4 
has focused on the effects of the hyperfine coupling 
and off-diagonal dipolar terms resulting in corrections to 
the above Hamiltoniani 14 i 15 i 16 i 17 Yet a non-perturabative 
calculation beyond mean-field of several fundamental 
properties such as the specific heat and linear suscep- 
tibility is lacking even for the first-order model described 
by Eq. {TJ. The goal of the present work is to numeri- 
cally investigate the above model and determine to what 
extent it can be used to interpret the experimental re- 
sults. In particular we calculate the static susceptibility 
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and specific heat and compare our result to recent experi- 
ments. We also calculate the ground state magnetization 
as a function of temperature in order to get an inde- 
pendent estimate of the critical dilution, x c , where the 
magnetization vanishes. Finally, we reexamine the low- 
temperature disordered phase and search for evidence of 
a stable spin-glass phase. 

II. METHOD 

We have used a single spin flip Monte Carlo method 
and applied periodic boundary conditions. To handle the 
long-range nature of the interaction we have used the 
Ewald summation method^ as explained in an earlier 
study£ To overcome energy barriers in the glassy phase 
and reach lower temperatures than in previous work we 
have used the replica exchange Monte Carlo method^ 
The method involves simulating an ensemble of systems 
at suitably chosen temperatures T il and the algorithm 
has two main phases. In the first phase each replica is 
independently evolved in (Monte Carlo) time using the 
single spin Metropolis algorithm. In the second phase 
attempts are made to exchange the replicas at adjacent 
temperatures Tj and Tj+x- A full Monte Carlo step con- 
sists of one attempted spin flip per spin (on average) 
followed by ten attempts to exchange neighboring repli- 
cas. For the simulation to converge at low temperatures 
it is important that the swap rate of the replicas is not 
too low. Theory and empirical studie a 20 i 21 have shown 
the optimal rate to be around 20% and these simulations 
were carried out with swap rates > 20% . 

In the present study we have calculated the specific 
heat, 

c =k^^ H2 )-m 2 )' ( 2 ) 

and the magnetic susceptibility 

X = - ((M 2 ) - <Af) 2 )) , (3) 

where the magnetization M = Yli=i °f ■ 

In order to study the disordered phase we have 
calculated the Edwards- Anderson overlap between two 
replicas* 2 ^ 

1 \ z z 

i 

and the corresponding Binder ratio 
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The spin glass susceptibility is defined as xsg = (q 2 )/T 2 . 
In addition to the thermal average we have calculated an 
average of 400-600 quenched disorder realizations. 



III. RESULTS 
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FIG. 1: Specific heat as a function of temperature. Experi- 
mental data from Ref. [1(3 (open symbols) for dilution x=0.08, 
x=0.045 and x=0.018 (top to bottom). Monte Carlo results 
(solid lines) for dilution x=0.165, 0.12, 0.08, 0.045 and 0.018 
(top to bottom). To display the limited finite-size effects re- 
sults for 8 3 and 10 3 (upper curve) unit cells are shown for 
each dilution. 

First we will compare our calculation of the specific 
heat to experimental data. In Fig. [1] we show our results 
for the specific heat as a function of temperature. In 
the same figure we also show recent experimental data. 10 
There is qualitative agreement and both sets of curves 
indicate that the specific heat grows with decreasing di- 
lution. Both sets of curves exhibit a maximum for some 
intermediate temperature, but the experimental peak po- 
sition is roughly independent of the dilution, while the 
calculated peak position scales linearly with x, which can 
be seen if Fig. [5] In an Ising spin glass the specific heat 
exhibits a broad maximum in the vicinity of the tran- 
sition temperature,— and as the mean-field transition 
temperature scales linearly with dilution, we might ex- 
pect the scaling observed in our calculation. However, 
the present experimental data are for high dilution, and 
experimentally there is no spin glass transition in the 
limit of extreme dilution. It would therefore be of inter- 
est to measure the specific heat for less dilute systems, 
to see whether the expected linear increase in the peak 
position of the susceptibility is recovered in this limit. 
Furthermore, the total specific heat is dominated by the 
contribution from the nuclear spinsi - and the necessary 
subtraction required to obtain the experimental data is 
sensitive to the form of the subtracted single-ion contri- 
bution. Any uncertainty in the noninteracting specific 
heat could cause a big change in the resulting electronic 
specific heat. 

Next we will analyze our results for the linear suscep- 
tibility. The inverse susceptibility is plotted in Fig. [3] 
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FIG. 2: Curie- Weiss temperature (Monte Carlo), mean- 
field critical temperature, maximum in specific heat (Monte 
Carlo), and experimental Curie- Weiss temperature from 
Ref. 2 (top to bottom) as a function of dilution. 




FIG. 4: Susceptibility per spin as a function of temperature 
for x=0.167, 0.12, 0.08 and 0.045 (top to bottom, solid lines). 
The dashed lines have slope -1 and -0.75 (lower line). 
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FIG. 3: Inverse susceptibility per spin as a function of tem- 
perature for x=0.045, 0.08, 0.12 and 0.167 (top to bottom). 
The displayed data has converged in system size. 



In good qualitative agreement with early measurements 
of the susceptibility^ we see Curie- Weiss behavior, \ ~ 
1/(T — T cw ), at higher temperatures, and deviations at 
low temperatures. As the dilution is increased the sus- 
ceptibility approaches the free spin limit x ~ T as 
expected. Extrapolating the Curie scaling to the inter- 
cept gives us the Curie- Weiss temperature, Tcw, which 
is positive, in accordance with the ferromagnetic correla- 
tions in LiHoF4. As can be seen in Fig. [5] we find that 
Tcw scales linearly in x, and we get Tcw = 0.08 and 
0.30 for x=0.045 and 0.167 respectively. Experiments 
reported in Ref. [H found that Tcw = 0.05 and 0.16 for 
x=0.045 and 0.167 respectively. From Fig. [5] we see that 



while our calculated Curie- Weiss temperature is higher 
than the mean-field critical temperature, the experimen- 
tal results are lower. We therefore reach qualitative, but 
not quantitative, agreement with this set of experiments. 
However, there is no experimental consensus on the func- 
tional form of the susceptibility and a later set of mea- 
surements by the same group report a scaling of the form 
X ~ T~ 0,75 for x=0.045. s In order to further analyze the 
functional form we plot our results for the susceptibility 
in a log-log plot in Fig. [4] From the inserted straight 
lines we see the Curie scaling x ~ T Q with a = — 1 at 
higher temperatures. As the temperature is lowered the 
susceptibility diverges faster, with an exponent a < — 1, 
contrary to reported measurements a = — 0.75£ The ex- 
perimental data was explained by off-diagonal terms in 
the dipolar interaction that arise when the material is 
diluted. Our omission of these terms could explain the 
discrepancy, but the fact remains that our results agree 
quite well with the earlier measurements of the suscepti- 
bility. Furthermore, the reported scaling of x ~ T~ 75 
persists up to T = 2 K, and given that the average diag- 
onal local dipolar field is of the order 1.53 x 0.045 ~ 0.07 
K it is surprising that there are deviations from Curie 
scaling at such elevated temperatures. Finally, experi- 
mental data for the susceptibility has also been argued 
to be well modeled by an exponential, low-temperature 
forrn^ x = ex P( — T/T)). Plotting our results for x in a 
semilog plot does not result in a straight line over any 
significant temperature interval. 

In order to compare the various results for the static 
susceptibility for the extreme dilution x = 0.045 we dis- 
play all the measurements in Fig. We have shifted the 
curves vertically to displaythe functional form better. 
We see that data from Ref . @ follows the form x ~ T -0 ' 75 
over the whole temperature range from 0.05 K to 2K. In 
the high temperature limit our calculation, as well as data 
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from Ref. and Ref. d tend to the Curie scaling \ ~ T~ 1 . 
At low temperature our calculation and data from Ref. 
diverge faster than T _1 while data from Ref. H grows sig- 
nificantly more slowly. Due to this discrepancy between 
different experiments it is difficult to determine how well 
the classical dipolar Ising model reflects the magnetic be- 
havior of LiHo 2; Yi_ 2 ;F4 in the high-dilution limit. More 
measurements that could explain the above experimen- 
tal differences would be necessary in order to draw more 
definite conclusions. 
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FIG. 5: Susceptibility at x=0.045. From top to bottom the 
open symbols are experimental data from Ref . |2L Monte Carlo 
results, experimental data from Ref. 9 and Ref. The dashed 
lines have slope -1 and -0.75 (lower curve). The curves have 
been separated vertically. 

Next we consider the magnetization of the dilute 
model. Using a parallel tempering method we are able 
to determine the magnetization curves to lower temper- 
ature than in a previous study— as can be seen in Fig. [6] 
Notice that, for a given dilution, the magnetization in- 
creases with system size for low dilution (x — 0.375), 
while it decreases in the more dilute systems [x — 0.25 
and x — 0.30). Extrapolating to the ground state we 
obtain the ground state magnetization curve in Fig. [7] 
Finite size effects and statistical errors prevent us from 
a very exact determination, but the curve indicates that 
the the critical concentration x c , where the magnetiza- 
tion vanishes, is about x c = 0.24 ± 0.03, which is a bit 
higher than the value x c = 0.21 ± 0.02 reported in a 
previous calculation 5 -. We are not aware of any precise 
experimental determination of x c , but our result is con- 
sistent with experiments that report a spin-glass phase at 
x = 0.167, but a ferromagnetic state at x = 0.3~ Our re- 
sults also compare well with a previous zero-temperature 
Monte Carlo study of Ising dipoles on a diluted BCC 
lattice,— where it was found that x c = 0.3 ± 0.1. 

Finally we consider the disordered phase for x < x c . 
We calculate the Binder ratio for the spin overlap, g q , in 
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FIG. 6: Magnetization per spin as a function of temperature 
for dilution x=0. 25, 0.3, 0.33, 0.375, 0.5, 0.625, 0.75 and 1 (left 
to right). At low dilution the result is for 10 3 unit cells, while 
for x=0.3, 0.33 and 0.375 the system sizes are 6 3 (dashed line), 
8 3 (dotted line) and 10 3 unit cells. For the highest dilution 
(x=0.25) the system sizes are 12 3 (dashed line), 14 3 (dotted 
line) and 16 3 unit cells. 
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FIG. 7: Ground state magnetization per spin as a function of 
dilution. 



the disordered phase. If there is a stable glass phase the 
curves for different system sizes are expected to cross at 
the freezing temperature. In a previous study 5 - no cross- 
ing was found, indicating that there is no freezing of the 
spin glass. Here we have repeated the calculation using 
the parallel tempering method in order to obtain more 
reliable data in the highly disordered phase. The result 
is shown in Fig. and as can be seen we still detect no 
crossing for x = 0.167. In order to analyze the nature 
of the disordered phase further we also consider the spin 
glass susceptibility xsg- In a study of the Heisenberg 
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FIG. 8: Binder ratio for the spin overlap at dilution x=0.167. 
System sizes are 10 3 , 12 3 and 14 3 unit cells (top to bottom). 




FIG. 9: Inverse spin glass susceptibility at dilution x = 0.167. 
The displayed data has converged in system size. 



spin glass it was argued that the Binder ratio of the spin 
overlap may not intersect at the freezing temperature for 
all boundary conditions^ The study suggests that the 
divergence of the spin glass susceptibility may be a better 
indicator of the freezing transition. Since many proper- 
ties of the long-range dipolar model are quite sensitive 
to the choice of boundary conditions, we therefore show 
results for the inverse spin glass susceptibility in Fig. [5] 
The finite size effects are very small and the spin glass 
susceptibility does not diverge at a finite temperature. 
Since there is quite convincing experimental evidence for 



a spin-glass transition^ at x = 0.167 the results are puz- 
zling and either there are some aspects of the simulations 
of the glassy dipolar phase that differ from the short- 
range Ising spin glass, or the neglected off-diagonal terms 
in the Hamiltonian are necessary to stabilize the glassy 
phase observed in LiHoo.i67Yo.833F4- 



IV. DISCUSSION AND CONCLUSION 

In this work we have made direct comparisons between 
calculations done on a first-order effective Ising dipolar 
model and experimental data obtained for LiH0.jY1_.rF4. 
We have focused on the static susceptibility and specific 
heat. Probably due to the slow dynamics in the highly 
disordered phase, different sets of experiments do not 
agree very well. Nevertheless, our calculation agrees well 
with static susceptibility measurements of Ref. |2| in the 
highly disordered regime. Obtaining non-perturbative re- 
sults beyond mean-field for the first-order classical dipo- 
lar model is an essential step on the way to understanding 
the physical properties of LiHo^Yi^Fzi. Three partic- 
ularly puzzling experimental results that cannot be ex- 
plained by our calculation on the classical dipolar model 
are the susceptibility scaling \ ~ T -0 75 for x — 0.045, 
a specific heat maximum that is independent of dilution, 
and the finite temperature spin-glass transition. 

The difference between our calculation and the experi- 
mental results may be explained by quantum mechanical 
terms that are not included in our classical model. We 
have ignored the hyperfine coupling between nuclear and 
electronic spins. In the low-temperature limit this is gen- 
erally important, particularly in the presence of a exter- 
nal transverse field. In this study we do not consider an 
applied magnetic field, and the hyperfine coupling is only 
expected to re-normalize the inter-spin coupling. 16 Ne- 
glecting the hyperfine coupling should therefore not lead 
to quantitatively new results, but may explain some of 
the quantitative differences between our calculations and 
the experiments. We have also ignored off-diagonal terms 
in the dipolar Hamiltonian, which result in an effective 
random transverse fiekh 16 ' 17 ! 25 Including these terms in 
the calculation would be an important next step in inter- 
preting the measurements on LiHo a; Yi_ a ;F4. 
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